In [1]:
from sympy import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
x1, x2, x3 = symbols("x_1 x_2 x_3")
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3")
R, L, ga, gv = symbols("R L g_a g_v")
init_printing()
In [2]:
a1 = pi / 2 + (L / 2 - alpha1)/R
x = R * cos(a1)
y = alpha2
z = R * sin(a1)
r = x*N.i + y*N.j + z*N.k
In [3]:
r
Out[3]:
In [4]:
r1 = r.diff(alpha1)
r2 = r.diff(alpha2)
k1 = trigsimp(r1.magnitude())
k2 = trigsimp(r2.magnitude())
r1 = r1/k1
r2 = r2/k2
In [5]:
r1
Out[5]:
In [6]:
r2
Out[6]:
In [7]:
n = r1.cross(r2)
n = trigsimp(n.normalize())
n
Out[7]:
In [8]:
dr1=r1.diff(alpha1)
k1 = trigsimp(r1.cross(dr1).magnitude()/k1**3)
k1
Out[8]:
In [9]:
dr2=r2.diff(alpha2)
k2 = trigsimp(r2.cross(dr2).magnitude()/k2**3)
k2
Out[9]:
In [10]:
n.diff(alpha1)
Out[10]:
$ \frac { d\vec{n} } { d\alpha_1} = -\frac {1}{R} \vec{v} = -k \vec{v} $
In [11]:
r1.diff(alpha1)
Out[11]:
$ \frac { d\vec{v} } { d\alpha_1} = \frac {1}{R} \vec{n} = k \vec{n} $
$ \vec{u} = u_v \vec{v} + u_n\vec{n} $
$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u_v\vec{v}) } { d\alpha_1} + \frac { d(u_n\vec{n}) } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_1} + \frac { du_v } { d\alpha_1} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} - u_n k \vec{v} + \frac { du_v } { d\alpha_1} \vec{v} + u_v k \vec{n}$
Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du_v } { d\alpha_1} - u_n k \right) \vec{v} + \left( \frac { du_n } { d\alpha_1} + u_v k \right) \vec{n}$
$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u_n\vec{n}) } { d\alpha_2} + \frac { d(u_v\vec{v}) } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_2} + \frac { du_v } { d\alpha_2} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + \frac { du_v } { d\alpha_2} \vec{v} $
In [12]:
R_alpha=r+alpha3*n
R_alpha
Out[12]:
In [13]:
R1=R_alpha.diff(alpha1)
R2=R_alpha.diff(alpha2)
R3=R_alpha.diff(alpha3)
In [14]:
trigsimp(R1)
Out[14]:
In [15]:
R2
Out[15]:
In [16]:
R3
Out[16]:
In [17]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import cm
x = R * cos(a1)
z = R * sin(a1)
alpha1_x = lambdify([R, L, alpha1], x, "numpy")
alpha3_z = lambdify([R, L, alpha1], z, "numpy")
R_num = 1/0.8
L_num = 2
x1_start = 0
x1_end = L_num
x3_start = -0.05
x3_end = 0.05
plot_x1_elements = 100
dx1 = (x1_end - x1_start) / plot_x1_elements
X_init = []
Y_init = []
x2 = 0
x3 = 0
for i in range(plot_x1_elements + 1):
x1 = x1_start + i * dx1
x=alpha1_x(R_num, L_num, x1)
z=alpha3_z(R_num, L_num, x1)
X_init.append(x)
Y_init.append(z)
plt.plot(X_init, Y_init, "r", label="початкова конфігурація")
geometry_title = "K={}".format(1/R_num)
plot_title = r"Форма панелі $L={}, h={}$".format(x1_end - x1_start, x3_end - x3_start)
if (len(geometry_title) > 0):
plot_title = r"Форма панелі $L={}, h={}, {}$".format(x1_end - x1_start, x3_end - x3_start, geometry_title)
plt.title(plot_title)
plt.axes().set_aspect('equal', 'datalim')
plt.legend(loc='best')
plt.xlabel(r"$x_1$, м", fontsize=12)
plt.ylabel(r"$x_3$, м", fontsize=12)
plt.grid()
plt.show()
In [37]:
eps=trigsimp(R1.dot(R2.cross(R3)))
R_1=simplify(trigsimp(R2.cross(R3)/eps))
R_2=simplify(trigsimp(R3.cross(R1)/eps))
R_3=simplify(trigsimp(R1.cross(R2)/eps))
eps
Out[37]:
In [19]:
R_1
Out[19]:
In [20]:
R_3
Out[20]:
$ A = \left( \begin{array}{ccc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \end{array} \right)$
$ \left[ \begin{array}{ccc} \vec{R}_1 & \vec{R}_2 & \vec{R}_3 \end{array} \right] = \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] \cdot \left( \begin{array}{ccc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \end{array} \right) = \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] \cdot A$
$ \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] =\left[ \begin{array}{ccc} \vec{R}_1 & \vec{R}_2 & \vec{R}_3 \end{array} \right] \cdot A^{-1}$
In [21]:
dx1da1=R1.dot(N.i)
dx1da2=R2.dot(N.i)
dx1da3=R3.dot(N.i)
dx2da1=R1.dot(N.j)
dx2da2=R2.dot(N.j)
dx2da3=R3.dot(N.j)
dx3da1=R1.dot(N.k)
dx3da2=R2.dot(N.k)
dx3da3=R3.dot(N.k)
A=Matrix([[dx1da1, dx1da2, dx1da3], [dx2da1, dx2da2, dx2da3], [dx3da1, dx3da2, dx3da3]])
simplify(A)
Out[21]:
In [22]:
A_inv = trigsimp(A**-1)
simplify(trigsimp(A_inv))
Out[22]:
In [23]:
trigsimp(A.det())
Out[23]:
${\displaystyle \hat{G}=\sum_{i,j} g^{ij}\vec{R}_i\vec{R}_j}$
In [24]:
g11=R_1.dot(R_1)
g12=R_1.dot(R_2)
g13=R_1.dot(R_3)
g21=R_2.dot(R_1)
g22=R_2.dot(R_2)
g23=R_2.dot(R_3)
g31=R_3.dot(R_1)
g32=R_3.dot(R_2)
g33=R_3.dot(R_3)
G=Matrix([[g11, g12, g13],[g21, g22, g23], [g31, g32, g33]])
G=trigsimp(G)
G
Out[24]:
${\displaystyle \hat{G}=\sum_{i,j} g_{ij}\vec{R}^i\vec{R}^j}$
In [25]:
g_11=R1.dot(R1)
g_12=R1.dot(R2)
g_13=R1.dot(R3)
g_21=R2.dot(R1)
g_22=R2.dot(R2)
g_23=R2.dot(R3)
g_31=R3.dot(R1)
g_32=R3.dot(R2)
g_33=R3.dot(R3)
G_con=Matrix([[g_11, g_12, g_13],[g_21, g_22, g_23], [g_31, g_32, g_33]])
G_con=trigsimp(G_con)
G_con
Out[25]:
In [26]:
G_inv = G**-1
G_inv
Out[26]:
In [27]:
dR1dalpha1 = trigsimp(R1.diff(alpha1))
dR1dalpha1
Out[27]:
$ \frac { d\vec{R_1} } { d\alpha_1} = -\frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \vec{R_3} $
In [28]:
dR1dalpha2 = trigsimp(R1.diff(alpha2))
dR1dalpha2
Out[28]:
In [29]:
dR1dalpha3 = trigsimp(R1.diff(alpha3))
dR1dalpha3
Out[29]:
$ \frac { d\vec{R_1} } { d\alpha_3} = \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \vec{R_1} $
In [30]:
dR2dalpha1 = trigsimp(R2.diff(alpha1))
dR2dalpha1
Out[30]:
In [31]:
dR2dalpha2 = trigsimp(R2.diff(alpha2))
dR2dalpha2
Out[31]:
In [32]:
dR2dalpha3 = trigsimp(R2.diff(alpha3))
dR2dalpha3
Out[32]:
In [33]:
dR3dalpha1 = trigsimp(R3.diff(alpha1))
dR3dalpha1
Out[33]:
$ \frac { d\vec{R_3} } { d\alpha_1} = \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \vec{R_1} $
In [34]:
dR3dalpha2 = trigsimp(R3.diff(alpha2))
dR3dalpha2
Out[34]:
In [35]:
dR3dalpha3 = trigsimp(R3.diff(alpha3))
dR3dalpha3
Out[35]:
$ \frac { d\vec{R_3} } { d\alpha_3} = \vec{0} $
$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u^1\vec{R_1}) } { d\alpha_1} + \frac { d(u^2\vec{R_2}) } { d\alpha_1}+ \frac { d(u^3\vec{R_3}) } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_1} + \frac { du^3 } { d\alpha_1} \vec{R_3} + u^3 \frac { d\vec{R_3} } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \vec{R_3} + \frac { du^2 } { d\alpha_1} \vec{R_2}+ \frac { du^3 } { d\alpha_1} \vec{R_3} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \vec{R_1}$
Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du^1 } { d\alpha_1} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \right) \vec{R_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + \left( \frac { du^3 } { d\alpha_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \right) \vec{R_3}$
$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u^1\vec{R_1}) } { d\alpha_2} + \frac { d(u^2\vec{R_2}) } { d\alpha_2}+ \frac { d(u^3\vec{R_3}) } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} + \frac { du^3 } { d\alpha_2} \vec{R_3} $
Then $ \frac { d\vec{u} } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} + \frac { du^3 } { d\alpha_2} \vec{R_3} $
$ \frac { d\vec{u} } { d\alpha_3} = \frac { d(u^1\vec{R_1}) } { d\alpha_3} + \frac { d(u^2\vec{R_2}) } { d\alpha_3}+ \frac { d(u^3\vec{R_3}) } { d\alpha_3} = \frac { du^1 } { d\alpha_3} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_3} + \frac { du^2 } { d\alpha_3} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_3} + \frac { du^3 } { d\alpha_3} \vec{R_3} + u^3 \frac { d\vec{R_3} } { d\alpha_3} = \frac { du^1 } { d\alpha_3} \vec{R_1} + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \vec{R_1} + \frac { du^2 } { d\alpha_3} \vec{R_2}+ \frac { du^3 } { d\alpha_3} \vec{R_3} $
Then $ \frac { d\vec{u} } { d\alpha_3} = \left( \frac { du^1 } { d\alpha_3} + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \right) \vec{R_1} + \frac { du^2 } { d\alpha_3} \vec{R_2}+ \frac { du^3 } { d\alpha_3} \vec{R_3}$
$\nabla_1 u^1 = \frac { \partial u^1 } { \partial \alpha_1} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}$
$\nabla_1 u^2 = \frac { \partial u^2 } { \partial \alpha_1} $
$\nabla_1 u^3 = \frac { \partial u^3 } { \partial \alpha_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) $
$\nabla_2 u^1 = \frac { \partial u^1 } { \partial \alpha_2}$
$\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$
$\nabla_2 u^3 = \frac { \partial u^3 } { \partial \alpha_2}$
$\nabla_3 u^1 = \frac { \partial u^1 } { \partial \alpha_3} + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}$
$\nabla_3 u^2 = \frac { \partial u^2 } { \partial \alpha_3} $
$\nabla_3 u^3 = \frac { \partial u^3 } { \partial \alpha_3}$
$ \nabla \vec{u} = \left( \begin{array}{ccc} \nabla_1 u^1 & \nabla_1 u^2 & \nabla_1 u^3 \\ \nabla_2 u^1 & \nabla_2 u^2 & \nabla_2 u^3 \\ \nabla_3 u^1 & \nabla_3 u^2 & \nabla_3 u^3 \\ \end{array} \right)$
In [36]:
u1=Function('u^1')
u2=Function('u^2')
u3=Function('u^3')
q=Function('q') # q(alpha3) = 1+alpha3/R
K = Symbol('K') # K = 1/R
u1_nabla1 = u1(alpha1, alpha2, alpha3).diff(alpha1) + u3(alpha1, alpha2, alpha3) * K / q(alpha3)
u2_nabla1 = u2(alpha1, alpha2, alpha3).diff(alpha1)
u3_nabla1 = u3(alpha1, alpha2, alpha3).diff(alpha1) - u1(alpha1, alpha2, alpha3) * K * q(alpha3)
u1_nabla2 = u1(alpha1, alpha2, alpha3).diff(alpha2)
u2_nabla2 = u2(alpha1, alpha2, alpha3).diff(alpha2)
u3_nabla2 = u3(alpha1, alpha2, alpha3).diff(alpha2)
u1_nabla3 = u1(alpha1, alpha2, alpha3).diff(alpha3) + u1(alpha1, alpha2, alpha3) * K / q(alpha3)
u2_nabla3 = u2(alpha1, alpha2, alpha3).diff(alpha3)
u3_nabla3 = u3(alpha1, alpha2, alpha3).diff(alpha3)
# $\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$
grad_u = Matrix([[u1_nabla1, u2_nabla1, u3_nabla1],[u1_nabla2, u2_nabla2, u3_nabla2], [u1_nabla3, u2_nabla3, u3_nabla3]])
grad_u
Out[36]:
In [37]:
G_s = Matrix([[q(alpha3)**2, 0, 0],[0, 1, 0], [0, 0, 1]])
grad_u_down=grad_u*G_s
expand(simplify(grad_u_down))
Out[37]:
$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_3 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \\ \nabla_3 u_2 \\ \nabla_1 u_3 \\ \nabla_2 u_3 \\ \nabla_3 u_3 \\ \end{array}
\left( \begin{array}{c} \left( 1+\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_1} + u^3 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \\ \left( 1+\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_2} \\ \left( 1+\frac{\alpha_3}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_3} + u^1 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ \frac { \partial u^3 } { \partial \alpha_1} - u^1 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) $
$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_3 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \\ \nabla_3 u_2 \\ \nabla_1 u_3 \\ \nabla_2 u_3 \\ \nabla_3 u_3 \\ \end{array}
B \cdot \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ \frac { \partial u^1 } { \partial \alpha_3} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ u^3 \\ \frac { \partial u^3 } { \partial \alpha_1} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) $
In [38]:
B = zeros(9, 12)
B[0,1] = (1+alpha3/R)**2
B[0,8] = (1+alpha3/R)/R
B[1,2] = (1+alpha3/R)**2
B[2,0] = (1+alpha3/R)/R
B[2,3] = (1+alpha3/R)**2
B[3,5] = S(1)
B[4,6] = S(1)
B[5,7] = S(1)
B[6,9] = S(1)
B[6,0] = -(1+alpha3/R)/R
B[7,10] = S(1)
B[8,11] = S(1)
B
Out[38]:
In [39]:
B_con = zeros(9, 12)
B_con[0,1] = 1
B_con[0,8] = (1+alpha3/R)/R
B_con[1,2] = 1
B_con[2,0] = -1/(R*(1+alpha3/R))
B_con[2,3] = 1
B_con[3,5] = S(1)
B_con[4,6] = S(1)
B_con[5,7] = S(1)
B_con[6,0] = -1/(R*(1+alpha3/R))
B_con[6,9] = S(1)
B_con[7,10] = S(1)
B_con[8,11] = S(1)
B_con
Out[39]:
In [40]:
q=(1+alpha3/R)
g_=(q*q)
koef=g_.diff(alpha3)
In [41]:
u_down_to_up=eye(12)
u_down_to_up[0,0]=g_
u_down_to_up[1,1]=g_
u_down_to_up[2,2]=g_
u_down_to_up[3,3]=g_
u_down_to_up[3,0]=koef
u_down_to_up
Out[41]:
In [42]:
simplify(B_con*u_down_to_up)
Out[42]:
$ \left( \begin{array}{c} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{array}
E \cdot \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_3 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \\ \nabla_3 u_2 \\ \nabla_1 u_3 \\ \nabla_2 u_3 \\ \nabla_3 u_3 \\ \end{array} \right)$
In [43]:
E=zeros(6,9)
E[0,0]=1
E[1,4]=1
E[2,8]=1
E[3,1]=1
E[3,3]=1
E[4,2]=1
E[4,6]=1
E[5,5]=1
E[5,7]=1
E
Out[43]:
In [70]:
N = 3
grad_u = zeros(N,N)
for i in range(N):
for j in range(N):
grad_u[i,j] = Symbol(r'\nabla_{} u_{}'.format(i+1,j+1))
grad_u
metric_tensor_marix_symbol = MatrixSymbol('g', N, N)
metric_tensor_marix = Matrix(metric_tensor_marix_symbol)
metric_tensor_marix
Out[70]:
In [77]:
e_nl_matrix = S(1)/S(2)*(grad_u*metric_tensor_marix*grad_u.T)
e_nl_matrix
Out[77]:
In [82]:
e_nl=zeros(6, 1)
e_nl[0,0] = e_nl_matrix[0,0]
e_nl[1,0] = e_nl_matrix[1,1]
e_nl[2,0] = e_nl_matrix[2,2]
e_nl[3,0] = 2*e_nl_matrix[0,1]
e_nl[4,0] = 2*e_nl_matrix[0,2]
e_nl[5,0] = 2*e_nl_matrix[1,2]
e_nl
Out[82]:
In [83]:
grad_u_g_symbol = MatrixSymbol('a', N, N)
grad_u_g = Matrix(grad_u_g_symbol)
grad_u_g
Out[83]:
In [86]:
a_nl_matrix = grad_u_g*grad_u.T
In [87]:
a_nl=zeros(6, 1)
a_nl[0,0] = a_nl_matrix[0,0]
a_nl[1,0] = a_nl_matrix[1,1]
a_nl[2,0] = a_nl_matrix[2,2]
a_nl[3,0] = 2*a_nl_matrix[0,1]
a_nl[4,0] = 2*a_nl_matrix[0,2]
a_nl[5,0] = 2*a_nl_matrix[1,2]
a_nl
Out[87]:
In [94]:
N = 3
grad_u_v = zeros(N*N, 1)
print(grad_u_v.shape)
for i in range(N):
for j in range(N):
index = i*N+j
# print("i = {}, j = {}, index = {}".format(i,j,index))
grad_u_v[index,0] = grad_u[j,i]
grad_u_v
Out[94]:
In [99]:
E_NL = zeros(6,9)
E_NL[0,0] = grad_u_g[0,0]
E_NL[0,3] = grad_u_g[0,1]
E_NL[0,6] = grad_u_g[0,2]
E_NL[1,1] = grad_u_g[1,0]
E_NL[1,4] = grad_u_g[1,1]
E_NL[1,7] = grad_u_g[1,2]
E_NL[2,2] = grad_u_g[2,0]
E_NL[2,5] = grad_u_g[2,1]
E_NL[2,8] = grad_u_g[2,2]
E_NL[3,1] = 2*grad_u_g[0,0]
E_NL[3,4] = 2*grad_u_g[0,1]
E_NL[3,7] = 2*grad_u_g[0,2]
E_NL[4,2] = 2*grad_u_g[0,0]
E_NL[4,5] = 2*grad_u_g[0,1]
E_NL[4,8] = 2*grad_u_g[0,2]
E_NL[5,2] = 2*grad_u_g[1,0]
E_NL[5,5] = 2*grad_u_g[1,1]
E_NL[5,8] = 2*grad_u_g[1,2]
E_NL
Out[99]:
In [101]:
E_NL*grad_u_v
Out[101]:
In [103]:
a_values=S(1)/S(2)*(grad_u*metric_tensor_marix)
E_NL = zeros(6,9)
E_NL[0,0] = a_values[0,0]
E_NL[0,3] = a_values[0,1]
E_NL[0,6] = a_values[0,2]
E_NL[1,1] = a_values[1,0]
E_NL[1,4] = a_values[1,1]
E_NL[1,7] = a_values[1,2]
E_NL[2,2] = a_values[2,0]
E_NL[2,5] = a_values[2,1]
E_NL[2,8] = a_values[2,2]
E_NL[3,1] = 2*a_values[0,0]
E_NL[3,4] = 2*a_values[0,1]
E_NL[3,7] = 2*a_values[0,2]
E_NL[4,2] = 2*a_values[0,0]
E_NL[4,5] = 2*a_values[0,1]
E_NL[4,8] = 2*a_values[0,2]
E_NL[5,2] = 2*a_values[1,0]
E_NL[5,5] = 2*a_values[1,1]
E_NL[5,8] = 2*a_values[1,2]
E_NL
Out[103]:
In [98]:
Q=E*B
Q=simplify(Q)
Q
Out[98]:
$u^1 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u\left( \alpha_1 \right)+\alpha_3\gamma \left( \alpha_1 \right) $
$u^2 \left( \alpha_1, \alpha_2, \alpha_3 \right)=0 $
$u^3 \left( \alpha_1, \alpha_2, \alpha_3 \right)=w\left( \alpha_1 \right) $
$ \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ \frac { \partial u^1 } { \partial \alpha_3} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ u^3 \\ \frac { \partial u^3 } { \partial \alpha_1} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) = T \cdot \left( \begin{array}{c} u \\ \frac { \partial u } { \partial \alpha_1} \\ \gamma \\ \frac { \partial \gamma } { \partial \alpha_1} \\ w \\ \frac { \partial w } { \partial \alpha_1} \\ \end{array} \right) $
In [45]:
T=zeros(12,6)
T[0,0]=1
T[0,2]=alpha3
T[1,1]=1
T[1,3]=alpha3
T[3,2]=1
T[8,4]=1
T[9,5]=1
T
Out[45]:
In [46]:
Q=E*B*T
Q=simplify(Q)
Q
Out[46]:
In [47]:
from sympy import MutableDenseNDimArray
C_x = MutableDenseNDimArray.zeros(3, 3, 3, 3)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
el = Symbol(elem_index)
C_x[i,j,k,l] = el
C_x
Out[47]:
In [48]:
C_x_symmetry = MutableDenseNDimArray.zeros(3, 3, 3, 3)
def getCIndecies(index):
if (index == 0):
return 0, 0
elif (index == 1):
return 1, 1
elif (index == 2):
return 2, 2
elif (index == 3):
return 0, 1
elif (index == 4):
return 0, 2
elif (index == 5):
return 1, 2
for s in range(6):
for t in range(s, 6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
el = Symbol(elem_index)
C_x_symmetry[i,j,k,l] = el
C_x_symmetry[i,j,l,k] = el
C_x_symmetry[j,i,k,l] = el
C_x_symmetry[j,i,l,k] = el
C_x_symmetry[k,l,i,j] = el
C_x_symmetry[k,l,j,i] = el
C_x_symmetry[l,k,i,j] = el
C_x_symmetry[l,k,j,i] = el
C_x_symmetry
Out[48]:
In [49]:
C_isotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)
C_isotropic_matrix = zeros(6)
mu = Symbol('mu')
la = Symbol('lambda')
for s in range(6):
for t in range(s, 6):
if (s < 3 and t < 3):
if(t != s):
C_isotropic_matrix[s,t] = la
C_isotropic_matrix[t,s] = la
else:
C_isotropic_matrix[s,t] = 2*mu+la
C_isotropic_matrix[t,s] = 2*mu+la
elif (s == t):
C_isotropic_matrix[s,t] = mu
C_isotropic_matrix[t,s] = mu
for s in range(6):
for t in range(s, 6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
el = C_isotropic_matrix[s, t]
C_isotropic[i,j,k,l] = el
C_isotropic[i,j,l,k] = el
C_isotropic[j,i,k,l] = el
C_isotropic[j,i,l,k] = el
C_isotropic[k,l,i,j] = el
C_isotropic[k,l,j,i] = el
C_isotropic[l,k,i,j] = el
C_isotropic[l,k,j,i] = el
C_isotropic
Out[49]:
In [50]:
def getCalpha(C, A, q, p, s, t):
res = S(0)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
return simplify(trigsimp(res))
C_isotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
c = getCalpha(C_isotropic, A_inv, i, j, k, l)
C_isotropic_alpha[i,j,k,l] = c
C_isotropic_alpha[0,0,0,0]
Out[50]:
In [51]:
C_isotropic_matrix_alpha = zeros(6)
for s in range(6):
for t in range(6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
C_isotropic_matrix_alpha[s,t] = C_isotropic_alpha[i,j,k,l]
C_isotropic_matrix_alpha
Out[51]:
In [52]:
C_orthotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)
C_orthotropic_matrix = zeros(6)
for s in range(6):
for t in range(s, 6):
elem_index = 'C^{{{}{}}}'.format(s+1, t+1)
el = Symbol(elem_index)
if ((s < 3 and t < 3) or t == s):
C_orthotropic_matrix[s,t] = el
C_orthotropic_matrix[t,s] = el
for s in range(6):
for t in range(s, 6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
el = C_orthotropic_matrix[s, t]
C_orthotropic[i,j,k,l] = el
C_orthotropic[i,j,l,k] = el
C_orthotropic[j,i,k,l] = el
C_orthotropic[j,i,l,k] = el
C_orthotropic[k,l,i,j] = el
C_orthotropic[k,l,j,i] = el
C_orthotropic[l,k,i,j] = el
C_orthotropic[l,k,j,i] = el
C_orthotropic
Out[52]:
In [53]:
def getCalpha(C, A, q, p, s, t):
res = S(0)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
return simplify(trigsimp(res))
C_orthotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
c = getCalpha(C_orthotropic, A_inv, i, j, k, l)
C_orthotropic_alpha[i,j,k,l] = c
C_orthotropic_alpha[0,0,0,0]
Out[53]:
In [54]:
C_orthotropic_matrix_alpha = zeros(6)
for s in range(6):
for t in range(6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
C_orthotropic_matrix_alpha[s,t] = C_orthotropic_alpha[i,j,k,l]
C_orthotropic_matrix_alpha
Out[54]:
$u^1=\frac{u_{[1]}}{1+\frac{\alpha_3}{R}}$
$\frac{\partial u^1} {\partial \alpha_3}=\frac{1}{1+\frac{\alpha_3}{R}} \frac{\partial u_{[1]}} {\partial \alpha_3} + u_{[1]} \frac{\partial} {\partial \alpha_3} \left( \frac{1}{1+\frac{\alpha_3}{R}} \right) = =\frac{1}{1+\frac{\alpha_3}{R}} \frac{\partial u_{[1]}} {\partial \alpha_3} - u_{[1]} \frac{1}{R \left( 1+\frac{\alpha_3}{R} \right)^2} $
In [55]:
P=eye(12,12)
P[0,0]=1/(1+alpha3/R)
P[1,1]=1/(1+alpha3/R)
P[2,2]=1/(1+alpha3/R)
P[3,0]=-1/(R*(1+alpha3/R)**2)
P[3,3]=1/(1+alpha3/R)
P
Out[55]:
In [56]:
Def=simplify(E*B*P)
Def
Out[56]:
In [57]:
rows, cols = Def.shape
D_p=zeros(rows, cols)
q = 1+alpha3/R
for i in range(rows):
ratio = 1
if (i==0):
ratio = q*q
elif (i==3 or i == 4):
ratio = q
for j in range(cols):
D_p[i,j] = Def[i,j] / ratio
D_p = simplify(D_p)
D_p
Out[57]:
In [58]:
C_isotropic_alpha_p = MutableDenseNDimArray.zeros(3, 3, 3, 3)
q=1+alpha3/R
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
fact = 1
if (i==0):
fact = fact*q
if (j==0):
fact = fact*q
if (k==0):
fact = fact*q
if (l==0):
fact = fact*q
C_isotropic_alpha_p[i,j,k,l] = simplify(C_isotropic_alpha[i,j,k,l]*fact)
C_isotropic_matrix_alpha_p = zeros(6)
for s in range(6):
for t in range(6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
C_isotropic_matrix_alpha_p[s,t] = C_isotropic_alpha_p[i,j,k,l]
C_isotropic_matrix_alpha_p
Out[58]:
In [59]:
C_orthotropic_alpha_p = MutableDenseNDimArray.zeros(3, 3, 3, 3)
q=1+alpha3/R
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
fact = 1
if (i==0):
fact = fact*q
if (j==0):
fact = fact*q
if (k==0):
fact = fact*q
if (l==0):
fact = fact*q
C_orthotropic_alpha_p[i,j,k,l] = simplify(C_orthotropic_alpha[i,j,k,l]*fact)
C_orthotropic_matrix_alpha_p = zeros(6)
for s in range(6):
for t in range(6):
i,j = getCIndecies(s)
k,l = getCIndecies(t)
C_orthotropic_matrix_alpha_p[s,t] = C_orthotropic_alpha_p[i,j,k,l]
C_orthotropic_matrix_alpha_p
Out[59]:
In [60]:
D_p_T = D_p*T
K = Symbol('K')
D_p_T = D_p_T.subs(R, 1/K)
simplify(D_p_T)
Out[60]:
$A=\frac {\theta}{2} \left( R + h_2 \right)^2-\frac {\theta}{2} \left( R + h_1 \right)^2$
In [61]:
theta, h1, h2=symbols('theta h_1 h_2')
square_geom=theta/2*(R+h2)**2-theta/2*(R+h1)**2
expand(simplify(square_geom))
Out[61]:
${\displaystyle A=\int_{0}^{L}\int_{h_1}^{h_2} \left( 1+\frac{\alpha_3}{R} \right) d \alpha_1 d \alpha_3}, L=R \theta$
In [62]:
square_int=integrate(integrate(1+alpha3/R, (alpha3, h1, h2)), (alpha1, 0, theta*R))
expand(simplify(square_int))
Out[62]:
In [73]:
S = simplify(D_p.T*C_isotropic_matrix_alpha_p*D_p*(1+alpha3/R)**2)
S
Out[73]:
In [64]:
W = simplify(D_p_T.T*C_isotropic_matrix_alpha_p*D_p_T*(1+K*alpha3))
W
Out[64]:
In [65]:
h=Symbol('h')
E=Symbol('E')
v=Symbol('nu')
W_a3 = integrate(W, (alpha3, -h/2, h/2))
W_a3 = simplify(W_a3)
W_a3.subs(la, E*v/((1+v)*(1-2*v))).subs(mu, E/((1+v)*2))
Out[65]:
In [66]:
A_M = zeros(3)
A_M[0,0] = E*h/(1-v**2)
A_M[1,1] = 5*E*h/(12*(1+v))
A_M[2,2] = E*h**3/(12*(1-v**2))
Q_M = zeros(3,6)
Q_M[0,1] = 1
Q_M[0,4] = K
Q_M[1,0] = -K
Q_M[1,2] = 1
Q_M[1,5] = 1
Q_M[2,3] = 1
W_M=Q_M.T*A_M*Q_M
W_M
Out[66]:
In [67]:
rho=Symbol('rho')
B_h=zeros(3,12)
B_h[0,0]=1
B_h[1,4]=1
B_h[1,8]=1
M=simplify(rho*P.T*B_h.T*G_con*B_h*P)
M
Out[67]:
In [74]:
M_p = rho*B_h.T*B_h*(1+alpha3/R)**2
M_p
Out[74]:
In [76]:
mass_matrix_func = lambdify((rho, R, alpha3), M_p, "numpy")
mass_matrix_func(100,10,20)
Out[76]:
In [77]:
stiffness_matrix_func = lambdify([R, mu, la, alpha3], S, "numpy")
stiffness_matrix_func(100, 200, 300, 400)
Out[77]:
In [78]:
import fem.geometry as g
import fem.model as m
import fem.material as mat
import fem.solver as s
import fem.mesh as me
import plot
def generate_layers(thickness, layers_count, material):
layer_top = thickness / 2
layer_thickness = thickness / layers_count
layers = set()
for i in range(layers_count):
layer = m.Layer(layer_top - layer_thickness, layer_top, material, i)
layers.add(layer)
layer_top -= layer_thickness
return layers
def solve(width, curvature, thickness):
layers_count = 1
layers = generate_layers(thickness, layers_count, mat.IsotropicMaterial.steel())
mesh = me.Mesh.generate(width, layers, N, M, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
# geometry = g.CorrugatedCylindricalPlate(width, curvature, corrugation_amplitude, corrugation_frequency)
geometry = g.CylindricalPlate(width, curvature)
# geometry = g.Geometry()
model = m.Model(geometry, layers, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
return s.solve(model, mesh, stiffness_matrix, mass_matrix)
def stiffness_matrix(material, geometry, x1, x2, x3):
return stiffness_matrix_func(1/geometry.curvature, material.mu(), material.lam(), x3)
def mass_matrix(material, geometry, x1, x2, x3):
return mass_matrix_func(material.rho, 1/geometry.curvature, x3)
# r=2
# width = r*2*3.14
# curvature = 1/r
width = 2
curvature = 0.8
thickness = 0.05
N = 100
M = 10
results = solve(width, curvature, thickness)
results_index = 0
plot.plot_init_and_deformed_geometry(results[results_index], 0, width, -thickness / 2, thickness / 2, 0)
#plot.plot_init_geometry(results[results_index].geometry, 0, width, -thickness / 2, thickness / 2, 0)
# plot.plot_strain(results[results_index], 0, width, -thickness / 2, thickness / 2, 0)
to_print = 20
if (len(results) < to_print):
to_print = len(results)
for i in range(to_print):
print(results[i].freq)
In [ ]: